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ABSTRACT 

We present a new algorithm that can reconstruct the full distributions of metric com¬ 
ponents within the class of spherically symmetric dust universes that may include a 
cosmological constant. The algorithm is capable of confronting this class of solutions 
with arbitrary data and opens a new observational window to determine the value of 
the cosmological constant. In this work we use luminosity and age data to constrain 
the geometry of the universe up to a redshift of z = 1.75. We show that, although cur¬ 
rent data are perfectly compatible with homogeneous models of the universe, simple 
radially inhomogeneous void models that are sometimes used as alternative explana¬ 
tions for the apparent acceleration of the late time universe cannot yet be ruled out. 
In doing so we reconstruct the density of cold dark matter out to z = 1.75 and derive 
constraints on the metric components when the universe was 10.5 Gyr old within a 
comoving volume of approximately 1 Gpc^. 
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1 INTRODUCTION 

Assuming that general relativity correctly describes grav¬ 
ity on cosmological scales the problem of reconstructing the 
cosmological metric can be posed as a characteristic initial 
value problem (CIVP). Spherically symmetric dust mod¬ 
els that may include a cosmological constant (also known 
as A Lemaitre-Tolman-Bondi (ALTB) models) have, after 
all gauge freedoms have been fixed, two free functional de¬ 
grees of freedom and a free parameter viz. the cosmologi¬ 
cal constant A. Ideally the two functional degrees of free¬ 
dom should be fixed using data generated by astrophys- 
ical models that do not presuppose any specific cosmol¬ 
ogy. This is a major issue for the observational cosmology 
programme Ellis et al. (1985); Stoeger et al. (1992a,b,c,d, 
1994); Maartens et al. (1996); Stoeger et al. (1997); Araujo 
& Stoeger (1999); Araujo, M.E. and Stoeger, W.R. (2009); 
Araujo & Stoeger (2009, 2010, 2011) with research per¬ 
taining to observations in non-symmetric universes ongoing 
Hellaby (2013, 2001). This paper however is not about ob¬ 
servations in non-symmetric universes. Our aim instead is 
to specify an algorithm that can determine the full distribu¬ 
tion of metric components in a spherically symmetric uni¬ 
verse given observational data. As such we do not make any 
attempt to justify how the data are obtained. A detailed in¬ 
vestigation into the model dependence of currently available 
data, as well as data expected from future surveys, has been 
relegated to an accompanying paper that will be released at 
a later stage. Note also that the algorithm we propose tests 
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for radial inhomogeneities on large scales but does not ad¬ 
dress the averaging or backreaction problems (for reviews of 
these topics see Clarkson et al. (2011); Buchert & RAd’sAd’- 
nen (2012) and references therein). 

After verifying that the algorithm performs as expected on 
simulated data we investigate what constraints can be de¬ 
rived from currently available luminosity (from the Union 
2.1 compilation Suzuki et al. (2012)) and age (from Cosmic 
Chronometers Moresco et al. (2012)) data. Both data sets 
are reported as functions of redshift and we convert the 
luminosity data into angular diameter distance data while 
age data are used to get a measure of the longitudinal ex¬ 
pansion rate. We have also used the best fit Planck value 
of Ho = 67.3 zb 1.2 kms“^Mpc“^ Ade et al. (2014). Thus 
we effectively assume that model independent data for D{z) 
and H\\[z) are available without questioning the validity of 
such data in ALTB cosmologies. We will further go on to 
show how data for the energy density of cold dark matter 
p{z) and the redshift drift can be incorporated into 

the algorithm. Since only two free functions are required to 
set the initial conditions of the CIVP some care has to be 
taken to incorporate more data without over-specifying the 
model. Note that all the initial data except the value of the 
cosmological constant can in principle be set with any two 
of D[z), H\\[z) and p{z). In section 2.2 we derive an expres¬ 
sion for the redshift drift and in section 4, we also show how 
the value of A can be inferred from redshift drift data. It 
has long been known that redshift drift data will be an in¬ 
valuable observable in cosmology (see for example Sandage 
(1962); Loeb (1998); Corasaniti et al. (2007); Ceng et al. 
(2015); Balbi & Quercellini (2007); Uzan et al. (2008)) but 
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as far as we are aware this is a new result, at least for ALTB 
models. 

The numerical integration scheme we employ to solve the 
Cl VP requires continuous functions as input. Since obser¬ 
vational data are reported at discrete values of the redshift 
we need a way to smooth data. Gaussian process regres¬ 
sion (GPR) provides a convenient semi-parametric Bayesian 
methodology for this purpose. The main difficulty with using 
GPR, or in fact any non-parametric smoothing algorithm, is 
making sure that the reconstructed functions obey certain 
physical constraints. As we explain in section 3.2, there ex¬ 
ists a relatively simple way to enforce all physical constraints 
if we choose a particular set to sample from viz. H\\ [z)^ p[z) 
and A. Once these have been specified the model is com¬ 
pletely determined by solving the Einstein field equations in 
observational coordinates. A likelihood can then be assigned 
to any set of initial conditions by comparing the solution to 
the observational data. This provides a means of performing 
inference on H\\[z)^ p[z) and A. However we stress from the 
onset that the algorithm we present requires redshift drift 
data to perform inference on A. Without data the 

best we can do is marginalise over the value of A by specify¬ 
ing some suitable prior. Prior specification is an important 
aspect of this algorithm that will be discussed in section 4. 
By now a number of algorithms have been proposed to solve 
the cosmological inverse problem (see Walters V Hellaby 
(2012); Redlich et al. (2014); Sapone et al. (2014); Valken- 
burg et al. (2014) and references therein). This paper ex¬ 
tends the formalism in Bester et al. (2014). To keep it as 
concise as possible we have omitted some details. In partic¬ 
ular we refer the reader to that paper, and van der Walt V 
Bishop (2012), for further details regarding our formulation 
of the observational and GIVP formalism. 

The paper is structured as follows: the next section high¬ 
lights some key differences between the ALTB and FLRW 
models that we exploit to test for radial inhomogeneity. We 
then give a brief description of the integration scheme used 
to numerically solve the Einstein field equations in observa¬ 
tional coordinates. To put our results within the standard 
cosmological context we also compute the coordinate trans¬ 
formation that relates observational coordinates to the stan¬ 
dard comoving 1+3 cosmological coordinates. In section 3 
we give a brief overview of how we employ Gaussian pro¬ 
cess regression to smooth the data on the current past light- 
cone (PLGO) with particular emphasis on how we enforce 
physical constraints. In section 4 we present the algorithm 
that allows us to reconstruct the distributions of the metric 
functions inside the past lightcone (PLG). We then verify 
that the algorithm performs as expected on simulated data. 
Finally in section 5 we present constraints from currently 
available data. We conclude with some prospects for future 
research. 

For notational convenience we denote 1 z = u through¬ 
out. A subscript zero can refer to a quantity evaluated on 
either the PLGO eg. D{wo,v) = Do{v) or the current time 
slice eg. H±{to,r) = H±o{r), the meaning should be unam¬ 
biguous from the context. The lower case Latin letter / is 
used to refer to probability distribution functions. Note in 
this paper an overdot refers to a partial derivative w.r.t. the 
coordinate w and a prime refers to a partial derivative w.r.t. 
the coordinate v and we work in units where c = G = 1 
throughout. 


2 FRAMEWORK 
2.1 The model 


The ALTB cosmological model describes a spherically sym¬ 
metric dust distribution Humphreys et al. (2012) in a uni¬ 
verse that may contain a cosmological constant Valken- 
burg (2012). In spherical synchronous comoving coordinates 
= [t,r, ^,0] (hereafter comoving coordinates or GG) the 
ALTB metric can be written as 


ds^ = —dt^ + X^{t, r)dr^ + R^(t, r)dQ?^ (1) 

where dO^ is the usual solid angle on the sphere, r is a co¬ 
moving radial coordinate and t measures proper time in the 
frame of the observer. One of the first integrals of the EFE’s 
allows us to relate X to R via X — g(r)drR where g(r) is an 
integration function that can be related to the space-time 
curvature. The EFE’s can further be used to write down the 
analogue of the Friedmann equation in ALTB models as 


(R) = 




-2 

flKO + f^AO, 


( 2 ) 


where H± = is the transverse (or perpendicular to the 
line of sight) expansion rate and the dimensionless density 
functions now depend or r. Explicitly they are defined by 




np 

3^’ 




A 

Ml 


and Qk = 1 — flm 


flA, (3) 


where p = p(t, r) is the energy density of cold dark matter 
and A is the cosmological constant. The expression (2) can 
be used to get the age of the universe as 


t - tB{r) 


1 dR^ 

iMJo rMKM)- 


(4) 


where tB{r) is a function of integration (the bang time) 
which can be scaled such that ts(0) = 0. In particular we 
get the current age of the universe at our space-time location 
by evaluating (4) with R = Ro along the central worldline 
of the observer (henceforth C) This fact will be used to 
set the time grid for the numerical integration scheme. 
Gomoving coordinates are probably the most intuitive co¬ 
ordinate system for cosmology. However they are not the 
most natural coordinate system when it comes to inter¬ 
preting data. For that we need observational coordinates 
= [w,v,0,(l)] (schematically illustrated in Figure 1). The 
metric in these coordinates can be written as 


ds^ = —A{w, v)dw^ + 2dwdv + D‘^{w, v)dQ? ^ (5) 


where re is a coordinate such that w = const, defines the 
PLG. Gauge fixing w to measure proper time along the cen¬ 
tral worldline of the observer then implies that A(w, 0) = 1. 
The coordinate v measures the null affine distance down 
the PLG and, with this form of the metric, is necessarily 
non-comoving. The coordinates 0 and 0 are spherical co¬ 
ordinates and are the same in the two coordinate systems. 
These definitions ensure that there is no gauge freedom left 
in these coordinates. This is important because it means 
we can specify any ALTB model by fixing two functional 
degrees of freedom and the parameter A. However it can 


^ Evaluated in practise by transforming this elliptic integral into 
one of the Carlson symmetric forms 
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Figure 1. A schematic illustration of the observational coordi¬ 
nate system. 


be difficult to interpret a cosmological model in terms of 
these coordinates mainly because they entangle spatial and 
temporal evolution. Therefore, to put our results in a more 
recognisable form, we compute the coordinate transforma¬ 
tion that sends observational to comoving coordinates. To 
do so requires specifying a gauge for r. We have chosen the 
gauge r = uD on the PLCO since it recovers the most well 
known form of the comoving FLRW metric viz. 

ds^ = —dt^ + a{t)^ (^ 1 . ( 6 ) 

ff the universe is truly FLRW then the comoving coordinates 
we transform to are exactly the same as those in the above 
metric, fn particular note that for a flat FLRW universe 
this gauge implies X(t^r) — a(t) and R{t,r) = a{t)r. Thus 
on any constant time slice the metric function in front of 
dr^ should be constant and the one in front of dQ^ should 
be a linear function of the comoving distance r. To further 
disentangle these two models we also identify functions that 
are a priori different in ALTB and FLRW universes. 


2.2 Discriminating between models 

A number of consistency relations have been proposed to 
test the assumption of homogeneity Maartens (2011). Here 
we review the fundamental differences between FLRW and 
ALTB models and show how they can be used to test for 
radial inhomogeneity. We will assume that D(z) and 
as well as all their derivatives w.r.t 2 : are known on the PLCO. 
We further assume that the metric components of (5) are 
known inside the PLC. The values of these quantities are 
all readily available once the CIVP has been solved (see 
section 2.3). We proceed with a discussion about kinematics 
in spherically symmetric dust space-time. 

First of all using UaU^ = —1 we write the four velocity as 
(recalling that = (1 + z) = u) 

47/2 _ 1 

= (7) 

Next a first order Taylor expansion on z{wo + Sw^vq + Sv) 
allows us, upon taking the limit of infinitesimal translations, 
to write 

dz _ dz . dz. dv . . 

^ “ ^l^^=const. + ^l-=const.^- («) 


Recognising the partial derivatives ^|-y=const. — u — 
if L=const. ~ u^H\\ and using the metric (5) to express 
^ we arrive at the following expression for the redshift drift 


dz 

dw 


_ • I f 2 tT 

= -u iL|| 



(9) 


fn section 4 we show how this fact can be used to infer the 
value of A. 

Next we can decompose the ray 4-vectors into parts parallel 
and orthogonal to as 

= (—Ufe/c^) (?x“ + e“) = — 2 /(?x“ + e“), (10) 

where = 1 and CaU^ = 0. 


Here e“ is the direction of propagation of the ray as seen by 

i.e. 


e = 


1 + , 
= 1“'^’- —’O’O]- 


(—Ubk^) ^ 2u 

Using this information to compute Hobs = |^ + 
gives (note that in ALTB models Hobs = = ^^) 

rr • uA-\-uA' 

Hobs = u+^+ ^ . ( 12 ) 

Further, projecting the null geodesic equation along the di¬ 
rection of propagation of the ray gives an expression that 
relates v to z viz. 


e^HVbka = 0 , 


dz 2 ^ 


(13) 


which can be used to get the affine parameter as a function 
of redshift 


.(^) = r 

Jo 


dz 


(14) 


0 {l + zYH^fzY 

The expression for the expansion scalar 0 = VaU^ is 
^ . uA' u'A u' 2uD D' uD'A 

The transverse expansion rate is found by noting that in 
spherical symmetry 0 = iLy + 2H^. Thus we have that 


(15) 


rr i ^ D' UAD' 


(16) 


The limiting behaviour for these expressions as u ^ 0 can be 
used to deduce that H^{w^t)) = H\\{w,0) i.e. the expansion 
rates are equal along C. This fact allows us to compute H± 
without resorting to series expansions. More importantly it 
allows to get the current age of the universe by evaluating 
(4) without direct observations of H±. In the remainder of 
the paper we use H to refer to the parallel expansion rate iLy 
and explicitly write H± for the transverse expansion rate. 
Note that the matter shear cr, which is identically zero in 
FLRW models, is related to the difference between H and 
H±. We can use this information to formulate a test for 
radial inhomogeneity as follows 

= (17) 

From the way Ti is defined it is clear that it is a dimension¬ 
less quantity that should be zero in FLRW models, ff it was 
possible to disentangle H from H± observationally we could 
formulate a very direct test based on this quantity (for a 
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related test see Maartens (2011), Clarkson et al. (2008)). 
Also, as proposed in Clarkson et al. (2008), we can formulate 
a test of the CP that is in principle independent of matter 
content and the particular theory of gravity employed by 
investigating the consistency between distances and the ex¬ 
pansion rate. The two main geometric effects on distance 
measurements are: 


(i) curvature bends null geodesics, 

(ii) expansion changes radial distances. 


In FLRW these two effects are coupled via the relation (as¬ 
suming distance duality (11 = where dz, is the lumi¬ 

nosity distance) 


D = 


1 


; sin J dz* 


Ho 


uHq^—Qko V Jo H{z*) 

This gives Qko in terms of H{z) and D{z) as 


Qko = 


[H {uD^, + D)f 
[HouDf 


(18) 


(19) 


Since Qro is expected to be independent of z this expression 
can be derived to yield a quantity that should be zero in 
FLRW models viz. 


T2 = 1 + [u{DD,,, - D%) - D^] 

T uHH^zD [iiD^z T , (20) 


where we use id = iLy since it is the radial expansion that 
effects radial geodesics. Note that although this quantity is 
in principle independent of the matter content and theory of 
gravity employed, the quantity that we reconstruct is not. 
This is because we use one of the field equations to constrain 
D(z) and its derivatives. Accurately correlating and recon¬ 
structing derivatives of D(z) and H(z) in a non-parametric 
way is a near impossible task otherwise. 

In what follows we will be using the quantities Ti and T 2 to 
measure the degree of radial inhomogeneity allowed by the 
data. Confirming that either of these quantities are zero, or 
at least small enough to be consistent with FLRW + pertur¬ 
bations, validates the Copernican principle against the data 
employed. 


2.3 Characteristic initial value problem 

Here we simply give a superficial outline of the CIVP for¬ 
malism and refer the reader to Bester et al. (2014) and van 
der Walt & Bishop (2012) for further details regarding its 
implementation. The field equations are 


D" = 

-^KDp(uif, 


(21) 

D' = 

^ [1 - DD'A' - 2DD' - A{D'f 



-ADD" - ^KpD^ - AD^ 


(22) 

A" = 

n' A'n' 

(23) 


D{0) = A'{0) = D{0) = 0, A(0) = D'{0) = 1. (24) 

These can be considered as constraint equations that need 
to be solved on each PLC. Equation (21) in particular plays 
a very important role in setting up the data on the PLCO 
since it ensures that all physical constraints are satisfied. 


Assuming that we have the functional forms of H(z) and 
p{z), as well as the value of A, the ALTB solution on the 
PLCO can be found with Procedure 1. This is all that is 


Procedure 1 _ 

1 H(z) with (14) gives the v{z) relation which gives 
p{v) and since v(z) is invertible also 

ui{v) = 1 + z{v); 

2 Solve (21) with ui{v) and p{v) to find D{v); 

3 Use (4) to get to and write initial data to the grid 
on which the solution is desired; 

4 Solve the coupled system (22) and (23) for A{v) and 
D{v) and use (25) to evaluate iii. 


required to compute the potential function 4> of (50) and 
thus confront the current realisations of H(z),p(z) and A 
with data. 

To find the solution on the inside of the PLC we need a way 
to evolve the initial data to the next PLC. To do so we need 
to specify, in addition to the output from Procedure 1, the 
values of uo,ui and p. The value of uo can be found from 
the normalisations condition UaU^ = — 1 whereas iii and p 
follow from the conservations equations VhT^a — 0- 


ui = 




(25) 


Ui 

.D 


D 


p{u\ + (u-o) T 2A u\ + 2—ui + A('ni) 


+ 2-^(ito + Aiti)) + p (ito + Aui) 


(26) 


There is still one final piece of information required to solve 
this system. As mentioned in section 2.1 the coordinate v 
is necessarily non-comoving. This means that its maximum 
extent changes as we go from one PLC to the next, its maxi¬ 
mum extent being determined by the causal horizon or char¬ 
acteristic cut off line (henceforth W). In order to compute 
this characteristic cut off we compute the path of a null 
geodesic from the maximum value of v into the interior of 
the PLC (see the discussion in Bester et al. (2014) for further 
details about how this is implemented). In principle anything 
beyond W has never been in causal contact with the observ¬ 
able universe. In a numerical application however, since at 
least two grid points are required to take derivatives, there 
is always a tiny portion of the PLC that is beyond our reach 
i.e. the intersection of W with C. This limitation is unavoid¬ 
able but note that the grid can be chosen fine enough so that 
these two grid points lie well within the averaging scale and 
are therefore of no practical importance. The numerical error 
of the integration scheme has been set to 10“® throughout. 
Since the scheme is second order the affine parameter dis¬ 
tance between these two grid points is v ^ ^10“® Gpc = 1 
Mpc. 

Finally, given the output from Procedure 1 above, the full 
ALTB solution in observational coordinates can be found 
using Procedure 2. To find the solution in comoving coor¬ 
dinates we need to specify the coordinate transformation. 
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Procedure 2 _ 

1 Get uq from UaU^ — —1 and evolve the system to 
the next PLC by solving (25) and then (26); 

2 Solve (21), (22) and (23) on this PLC; 

3 Repeat until the domain of calculation is exhausted; 

4 Get the characteristic cut off line W. 


2.4 Coordinate Transformation 


Here we derive the transformation that links observational to 
comoving coordinates (schematically illustrated in Figure 2). 
li t — t(r) is the solution for t on the PLCO then choosing 
R(t{r),r) = D(wo,v) ensures that 0 and 0 have the same 
meaning in both coordinate systems. Thus we only need to 
consider the (ic, v) ^ (t, r) transformation. The CIVP for¬ 
malism solves for the background cosmological metric (5) 
in terms of observational coordinates . We would like to 
compare these solutions to known solutions of the metric (1) 
in terms of comoving coordinates Accordingly we need 
to find both the metric components of (1) as functions of 
and then explicitly solve for comoving coordinates in terms 
of The metric components follow from the transforma¬ 
tion law 


_ dx^ dx^ 


(27) 


Glearly we need expressions for these partial derivatives 
purely in terms of the observational metric and coordinates. 
Once we have the comoving metric in observational coor¬ 
dinates the geodesic equations can be solved on each PLG 
to explicitly find comoving coordinates in terms of observa¬ 
tional coordinates. 

Gauge fixing the coordinate w to measure proper time along 
the central worldline means that the partial derivatives in¬ 
volving time are straightforward. Using that dt — —Uadx^, 
dw = kadx^ and = v,aU^ = ^ gives (where we have 
used the expression for from (7)) 


dt 

dw 

dw 

dt 


Av? + 1 
2u 


u and 


dt 

dv Av? — 1 
dt 2u 


(28) 

(29) 


The partial derivatives involving r require a little more ef¬ 
fort. The transformation (27) as well as its inverse gives 
(where o=|^,6=|^,c=|j) and d = for notational 
simplicity) 


d^X‘^ = g^d?{ab + hD'f = u, 

bu — l-auA — = 0, 

2 2u 

1 , . d ^ 

cu H —duA -= 0, 

2 2u 

ac-\-hd = 1, 


(30) 

(31) 

(32) 

(33) 


where we have used X — g(r)drR — g{r) [aD + hDj . Since 
there are five unknowns in four equations some additional 
information is required to solve this system. This is provided 
by the fact that the partial derivatives in these transforma¬ 
tions commute. No new information can be derived from 
applying this criterion to a and 5, it simply recovers the 
momentum conservation equation (25). Applied to c and d 
however we get a partial differential equation (PDF) for d 



Figure 2. An illustration of the intersection of the two coordinate 
systems. It is only possible to reconstruct constant time surfaces 
that lie within the 2 sphere of intersection of the coordinate sys¬ 
tems. 


viz. c = d. This casts (32) into the following flux conserva¬ 
tive form 




(34) 


Since A and u are known on the grid all that is required 
to solve this PDE are initial conditions for d. These are 
provided by fixing the gauge freedom in r. It is clearest to 
proceed with an FLRW analogy. To that end recall that in 
FLRW universes the comoving distance is simply r = 
where ac is a constant equal to the value of the scale factor 
on C. Thus in FLRW the transformation is extremely simple, 
d is always given by 

dr u'D + uD' ^ 

— =d= -. (35) 

dv ttc 

Having fixed d we may use (30)-(33) to find the remaining 
unknowns viz. a,b,c and X. To find the transformation in 
general we should see (35) as a normalisation on the PLGO. 
The initial data for d then makes it possible to solve (34) 
and the remaining unknowns again follow from (30)-(33). 
Partial derivatives w.r.t. either r or t can now be expressed 
as derivatives w.r.t. w and v using straightforward tensor 
calculus. To get t(w,v) and r(w,v) we further need to solve 
the following geodesic equations on each PLG 


d^t 

dtX dt 2 _ 

(M 

1 

(36) 

dv‘^ 

1 

d^r 

dv^ 

drX Ur\ 
X [dvj 

dr dt 
dvdv' 

(37) 


These solutions allow us to associate the corresponding (t, r) 
pair to any {w,v) grid point. In what follows we are going 
to reconstruct a snapshot of the universe at a fixed time, 
t = t* say, in the past. Let us denote the hypersurface de¬ 
fined by t* = const, by St and the radial coordinate on this 
hypersurface by r*. Then our aim is to find the intersection 
of functions defined on the {w,v) grid with St. This can be 
achieved by finding the intersection of each PLG for which 
w > t* with St and using the coordinate transformation 
to associate these points of intersection, (ic*,u*) say, to the 
corresponding comoving coordinates. This is implemented 
using Procedure 3. 
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Procedure 3 _ 

1 Use (28) and (29) to get partial derivatives w.r.t t ; 

2 Solve (34) for d with initial data set by (35) and use 
(30)-(33) to find a,5, c and X; 

3 Solve (36) and (37) on each PLC to get t{w,v) and 
r{w,v) on the grid.; 

4 By construction w* consists of all the grid points for 
which w > t*. Thus for each ic* find u* such that 
t(u*) = t* and then interpolate r(v) to find 

r* — r{y*)] 

5 Having located (ic*,u*) simply interpolate the 
desired functions and report their values at the 
corresponding (t*, r*). 

3 SMOOTHING THE DATA 
3.1 Gaussian process regression 

Gaussian processes have recently become quite popular as a 
means to perform non-parametric regression on cosmologi¬ 
cal data sets (see for example Seikel et ah (2012); Shafieloo 
et ah (2012); Busti et ah (2014a,b)). As a result we will only 
give a superficial outline of the theory relevant to the cur¬ 
rent application. In its most basic form a GP is a collection 
of random variables. Any finite collection of these variables 
have a joint Gaussian distribution Rasmussen X Williams 
(2006). A Gaussian process can be completely characterised 
by specifying its mean m(x) and covariance function /c(x, x). 
The mean and covariance function of a real process /(x) are 
defined by 

m(x) = E[/(x)], (38) 

fc(x,x) = E[(/(x)-m(x))(/(x)-m(x))], (39) 

where E[x] denotes the expectation value of x with respect to 
a Gaussian distribution. This is conveniently abbreviated us¬ 
ing the notation /(x) ~ 0P(m(x),/c(x, x)). We will further 
abbreviate our notation by labelling the Gaussian processes 
for the different functions using subscripts. Thus QVh and 
QVp refer to the Gaussian processes for the expansion rate 
and density respectively. 

The Gaussian process regression (GPR) problem consists of 
making inferences about the relationship between the in¬ 
puts and the targets. The covariance matrix follows from 
evaluating the covariance function at the relevant points i.e. 
Kij — k{xi,Xj). Denoting the prior distribution over func¬ 
tions as /p, the joint distribution between fp and the data 
is given by 


y 

~v(o, 

■ k{x,x) + j: 

K{X,Xp) 

fp 

KiXp,X) 

K{Xp,Xp) 


where S is the covariance matrix the of the data, X is the 
vector of inputs and Xp the vector of targets. The posterior 
distribution is the distribution of fp restricted to be com¬ 
patible with the observations or, in probabilistic terms, the 
conditional distribution of fp given data and targets i.e. 

fp I VyWp ~-V(/p,cov(/p)) , where (41) 
fp = K{Xp,X)Kf^{y-m{X)), (42) 

cov(/p) = K{Xp,Xp)-KiXp,X)Kf^KiX,Xp).i43) 

Here fp = E[/p|A, y, Ap] is the posterior mean, cov(/p) the 
posterior covariance matrix and Ky = [K(X,X) + E]. The 


marginal log-likelihood associated with GPR is given by 

\o^{p{v\X,e)) = -\y'^K-\- i log \Ky\- I log(2^), (44) 

where ^ is a vector of hyperparameters. Eqns (41) - (43) are 
the key predictive equations for GPR. 

We employ the Mattern class of covariance functions with 
v — 5/2 throughout (note r — \x — x\) 

fc(a:,®) = (l + ^^ + ^) cT/exp (-^^) . (45) 

Mean functions are chosen using an iterative procedure that 
will be discussed in section 4.1. 

3.2 Enforcing physical constraints 

Substituting = vfH into the chain rule applied to D' 
shows that the angular diameter distance and the expansion 
rate are related by 

, dzdD 

^ (46) 

dv dz 

D” = . (47) 

Our reasons for showing these equations are twofold. Firstly 
they show that incorporating the interdependence between 
D(z) and H(z) while smoothing them separately is non¬ 
trivial. It is far easier to start with realisations of H(z) and 
p{z) and then to find D{z) as described in Procedure 2. 
We then infer D(z) from the available data as explained in 
Algorithm 1 below. Note that regularity at the vertex of 
the cone is enforced by solving (21) with the specified ini¬ 
tial conditions (24). Secondly these expressions can be used 
to minimise the number of additional numerical derivatives 
that need to be found in order to reconstruct T 2 of (20). In 
fact we only need to numerically differentiate H w.r.t. ; 2 , the 
other required derivatives (towards z) then follow from (46) 
and (47). 

The null energy condition can be enforced simply by reject¬ 
ing samples with p < 0. To avoid shell crossing singularities 
we also need to ensure that the density remains regular as 
drR{t,r) crosses zero Krasihski (2014). This can only be 
done post integration and once we have the form of the co¬ 
ordinate transformation. We explicitly check this condition 
and discard the sample if it does. However since we sample 
p(z) directly we found that this virtually never happens. 


4 THE ALGORITHM 

Let us again emphasise that a spherically symmetric dust 
universe has two free functions. In our algorithm we select 
p(z) and H(z) as the two free functions. Given a value of A 
we then solve the system of field equations by implementing 
Procedure 1 followed by Procedure 2. Since the entire solu¬ 
tion is known numerically on a grid, and we can convert any 
function to a function of redshift, we are able to assign a 
likelihood to the current sample of H(z ), p{z) and A by con¬ 
fronting this solution with the available data. Importantly 
we can use any of the available data sets to compute this like¬ 
lihood. However we should point out that using only a single 
data set will tend to favour models of the universe that only 
have a single free function. It is therefore imperative that we 
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use at least two data sets to perform inference. Further, as 
evident from (21), density data can be used to further con¬ 
strain D{z) and H{z) but it cannot tell us anything about 
A since A has not yet entered the picture. Also, from (22), 
(23) and (25), the computation that leads to the form of 
ui{z) involves all of D{z), H{z), p{z) and A. Redshift drift 
data will therefore have the most constraining power and is 
required to ultimately test the Copernican principle. In full 
awareness of this limitation the best we can currently do is 
to marginalise over the value of A by sampling it from some 
suitable prior. As already mentioned prior specification is 
an important aspect of the algorithm so we will discuss that 
before presenting the algorithm used to perform inference. 

4.1 Prior specification 

Depending on the data available we might wish to specify 
priors in a number of different ways. The efficiency of our al¬ 
gorithm is vastly improved by knowing approximately which 
functions to look for. Specifically we use data to inform our 
priors when such data are available. In the absence of data 
we use our best guess and fine tune it by looking at the 
diagnostics of the MCMC. Of course there are many other 
ways in which priors could be specified, there is no claim 
that our choices are optimal in any sense. Actually a num¬ 
ber of improvements, eg. sieve priors Cotter et al. (2013), 
are possible. These will be investigated in future research. 
Note that the priors we specify below will also be used as 
proposals for the MCMC described in section 4.2. 

4 . 1.1 Sampling K 

Prior specification for the parameter A should be simple 
when redshift drift data become available. Since A is only 
a parameter any ID distribution should suffice. It would be 
simplest to incorporate a Gaussian prior into our algorithm 
because then we could simply augment the vector x in (49) 
with A and use the MCMC diagnostics to pin down its mean 
value and fine tune its variance. 

Since we cannot infer the value of A without redshift drift 
data we treat it as nuisance parameter. Depending on the 
application we use two separate priors. Firstly we do not 
want the uncertainty in A to obscure our results on simu¬ 
lated data. We therefore choose an accurate Gaussian prior 
when our primary intention is to verify that the algorithm 
works as expected on simulated data. In particular we use 
the background ACDM value with 2% uncertainty. 

When working with real data however we want to derive 
constraints regardless of the value of A. To do this we sam¬ 
ple A from a uniform distribution A ~ W[0,Amacc] where 
Amacc = 3(i/^)^. Here {Hq) is the 2-cr upper bound on the 
posterior of the expansion rate at the origin. Thus hmax is 
the 2-cr upper bound on A when Qa = I. 

4 . 1.2 Sampling H{z) 

There are a number of things to consider when setting a prior 
over functions for which data are available, two of which 
stand out in particular. Firstly, in order to avoid model mis- 
specification as much as possible, we do not want to presup¬ 
pose a parametrisation that might bias the space of functions 


considered. Secondly, it is important to explore the widest 
possible space of functions compatible with the data. With 
this in mind we used the following iterative procedure to set 
the prior on H(z). 

We start by performing GPR on expansion rate data using 
a zero mean function. We then use GaPP Seikel et al. (2012) 
to train the hyperparameters of the Gaussian process. The 
optimised hyperparameters together with (42) give the pos¬ 
terior mean function. We then redo the GPR using this as 
the prior mean function and iterate until there is no signif¬ 
icant change in the prior and the posterior mean functions. 
The posterior mean function for H(z) resulting from the fi¬ 
nal iteration is denoted by H. 

In the Bayesian sense the widest space of functions compat¬ 
ible with the data actually results from marginalising over 
the hyperparameters. However we found that, especially for 
small data sets, better results were obtained on simulated 
data by using the optimised values resulting from the final 
iteration above. These values are therefore used to compute 
the mean function and covariance matrix with (42) and (43) 
respectively. We then use (41) to draw smooth function re¬ 
alisations and use these as input to the CIVP. 

4 . 1.3 Sampling p(z) 

Given density data, the procedure outlined above to sample 
H{z) could simply be adapted to set a prior over p{z). In the 
absence of data we have to be explicit about how the prior 
for p(z) is set. Firstly we need to specify a mean function. 
For this we have implemented another iterative procedure in 
which we initially set the mean function as the background 
AGDM value. We then go through a number of burnin peri¬ 
ods each time using the posterior median of p(z) as the prior 
mean function in the next. This is repeated until there is no 
significant change in posterior median of p{z). We have ver¬ 
ified that specifying the initial mean function differently, as 
that corresponding to a void LTB model for instance, does 
not lead to a significantly different posterior median for p{z). 
The posterior median of p(z) resulting from the hnal itera¬ 
tion is denoted p. 

Of course to implement this iterative procedure we need to 
model the prior covariance matrix of p{z). For this we use a 
Gaussian process prior and then condition on artificial ob¬ 
servations Zi,p, 5Pi, where p gets updated after every burnin 
period. Within the GPR framework such a covariance ma¬ 
trix is given by 

T.j, = Kpp + Kj,{K+ (48) 

where E is a matrix with variances along the diagonal. We 
have chosen to model these artificial observations so that 
their variance increases as a power law of the redshift i.e. 
5p^ = a‘f {1 Zi)^. Here a is a real number chosen in a very 
‘liberal’ manner so that the 6pi completely overshoot the 
expected variance of p{z). It is clear from Figure 6 that this 
is indeed the case. 

The hyperparameters for GVp are chosen in a somewhat ad 
hoc manner. For the length scale I we initially pick a value 
that leads to reasonable density profiles. The signal vari¬ 
ance cFf is then set to some small value and increased until 
we see that the prior distribution of p{z) sufficiently covers 
its posterior. Both the length scale and the signal variance 
are adjusted iteratively to control the acceptance rate of the 
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Figure 3. Observables on the PLCO from simulated data. Top Left: Posterior distribution of D(z) on the PLCO. The contours are of 
similar width as when smoothing the D(z) alone except at the far end of the data where they are drastically reduced. This is because 
concave realisations of D(z) are excluded on the basis that they lead to violations of the null energy condition. Top Right: Posterior 
distribution of H(z) on the PLCO. The contours are significantly narrower than when smoothing H(z) data alone, especially for small 2 ;. 
The reason for this is that, as can clearly be seen from (46) and (47), the D(z) data are able to significantly constrain H{z). Bottom Left: 
The posterior distribution of p{z) on the PLCO. Clearly the error bars far overestimate the posterior uncertainty in the reconstructed 
p{z). Bottom Right: The posterior distribution of to- Note how close the peak is to the background ACDM value indicated by the black 
cross in the figure. Also to is quite well constrained, virtually none of the models have to < t*. 


functional MCMC detailed in Algorithm 1. 

This procedure might appear obscure and ad-hoc. However 
we found our results to be quite robust against this prior. 
Slightly altering the values of I and a/ can effect on the over¬ 
all acceptance rate of the MCMC but it does significantly 
change the final distributions of quantities Ti and T 2 used 
to test the Copernican principle, only the time it takes for 
these distributions to converge. 

4.2 Inference 

Here we describe the method used to perform inference. An¬ 
gular diameter distance and expansion rate data only allow 
us to perform inference on H[z) and p{z) not on A. As such 
we construct the target vector 


where H is the posterior mean function oi QVh and p the 
posterior median function of QVp. The quantity we want to 
infer from the available data is the posterior oi x. To do 
so we implement a modified random walk method over the 
function space of a:. A detailed discussion of MCMC algo¬ 
rithms for functional spaces is given in Cotter et al. (2013). 
We will restrict the current discussion only to the relevant 
theory. Note that pi{x) is used to refer to the measure on x 
while /io (a:) is used to refer to the measure on the prior over 

X. 


One of the key difficulties with MCMCs over function spaces 
is that their acceptance rates are sensitive to spatial discriti- 
zation. As the grid is refined acceptance rates typically drop. 
This is obviously far from ideal for numerical solutions to dif¬ 
ferential equations. The key idea to overcome this difficulty 
is to use, as proposals for Metropolis-Hastings type methods, 
discritizations of certain differential equations which exactly 
preserve the Gaussian reference measure as the likelihood 
drops to zero. In other words, in absence of evidence to the 
contrary, the posterior will be identical to the prior. Note 
that the prior p.o{x) is now the joint normal distribution of 
the priors over H and p. We assume that H and p are inde¬ 
pendent in the prior so that the off diagonal block matrices 
in the joint covariance matrix are just zero matrices. This 
effectively means that we can get a prior sample of x by 
sampling H[z) and p[z) separately. 

Next we note that if either the prior piQ{x) or the posterior 
/i(x) can be chosen to be a dominating reference Gaussian 
measure, then Bayes’ formula is expressed with the corre¬ 
sponding Radon-Nikodym derivative as 

-^(x)ocL(x), where L(x) = exp(—4>(x)). (50) 

a/io 

Here L is the likelihood with corresponding real valued po¬ 
tential 4> chosen so that the Markov chain converges to 
the desired target. The current application dictates that we 
choose the prior to be the dominating reference Gaussian 
measure. 
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Figure 4. The quantities Ti and T 2 from simulated data. Top left: Posterior distribution of Ti on the PLCO marginalised over to- In 
order to compare successive realisations we plot Ti as a function of the normalised affine parameter i.e. —-—. Middle Left: Posterior 

'^max 

distribution of Ti on the PLCF similarly plotted against the normalised affine parameter. Bottom Left: The maximum extent of v on 
the PLCF. Note how close the peak falls to the background ACDM value indicated by the black cross in the figure. Top Right: Posterior 
distribution of T 2 on the PLCO marginalised over to - In order to compare successive realisations we plot T 2 as a function of the normalised 
redshift i.e. —^—. Middle Right: Posterior distribution of T 2 on the PLCF plotted against the normalised redshift. Bottom Right: The 

^max 

maximum extent of 2 : on the PLCF. Again the peak is very close to the background ACDM value indicated by the black cross in the 
figure. 


The preconditioned Crank-Nicolson (pCN) proposal of Cot¬ 
ter et al. (2013) takes the form 

y(k) = V(1 - I3^)x(k) + hS, with 5~/io(a;), (51) 

where 0 < /3 < 1 is a constant that can be adjusted to 
control the acceptance rate and we indicate the step in the 
chain by a subscript in braces. The acceptance probability 
is then defined as 


a{x, y) = min (1, exp ($(a:) - $(y))). 


(52) 


This method differs from the standard random walk method 
since the proposal is not centred but rather autoregressive 
with order one {AR(1) type). Importantly the method does 
not reject the proposal in the case where $ = 0 but rather 
accepts the move with probability one (which is why prior 
specification is so important). This is the key idea, the 
method exactly preserves the Gaussian reference measure 
because the proposal is reversible w.r.t. /xo- If we want to 
ensure that the MCMC rejects a sample, for example when 
certain physical constraints have been violated, then con¬ 
ceptually we should set $ = oo. 

Next we need to specify the potential function 4>. For sim¬ 
plicity we use the joint chi-square of all the available data 
for this purpose 

= (53) 


Here Xi is the chi-square value of data set i. For example if 
we have j observations of the function y then we use 

2 _ 1 (Vj -vY 

Xy n 2 ^ 


(54) 


'Vi 


We could use a potential function other than the chi- 
square if we wanted to. Actually there is considerable flex¬ 
ibility in the way 4> can be specified and this is important 
for the intended application of the algorithm. As pointed 
out in the introduction it is substantially more difficult to 
obtain data that are valid in inhomogeneous cosmologies. 
This limits the applicability of our algorithm since there 
are very few model independent data available. For exam¬ 
ple there is an inherent circularity when incorporating BAO 
Bassett & Hlozek (2010), weak gravitational lensing Bartel- 
mann h Schneider (2001), redshift-space distortion Perci- 
val et al. (2011), galaxy number count Gardner (1998) and 
galaxy cluster Allen et al. (2011) data since they often as¬ 
sume a background FLRW cosmology at some point or an¬ 
other. Moreover, because of the complexity behind many 
of these astrophysical modelling techniques, the exact point 
at which such assumptions enter the calculation can be ob¬ 
scured to non-experts in any particular field. However, if we 
are able to pinpoint exactly where these assumptions enter 
the picture, it might still be possible to incorporate these 
data sets in a meaningful way. This would be possible if we 
can identify and reverse the effect that these assumptions 
have on the potential function 4>. A similar idea is used to 
analyse GMB data in Vonlanthen et al. (2010) and Audren 
(2014). Probably the simplest conceivable example would be 
to marginalise over the particular choice of iLo when convert¬ 
ing the Union 2.1 distance modulus to angular diameter dis¬ 
tance data. However having admitted that current data are 
not able to place tight constraints on violations of the Goper- 
nican principle, we chose not to perform the marginalisation 
since that would only degrade the constraints further. Also 
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Figure 5. Posterior distributions on a constant time slice Et defined by t* = 10.5 Gyr for simulated data. We plot X, R and p as 
functions of the normalised comoving distance —-—. Top Left: As expected the metric function X is very nearly constant. Bottom Left: 

'^max 

The density on St plotted in units of the critical density pc of the background ACDM model used in the simulations. Clearly the most 
probable density profile is that of a homogeneous matter distribution. Interestingly it seems that void profiles are not ruled out by the 
data. Top Right: Again as expected the metric function R is very nearly a linear function of the comoving distance r. Bottom Right: The 
maximum extent of the radial coordinate r* on St. The peak of this distribution is very close to the background ACDM value indicated 
by the black cross in the figure. 


note that Hq is not the only nuisance parameter involved in 
analysing supernovae data Amanullah et al. (2010). 

It is unlikely that future surveys such as Euclid Amendola 
et al. (2013), DES Sanchez & the Des collaboration (2010) or 
the SKA Maartens et al. (2015) will produce data which are 
completely model independent. We might therefore have to 
content ourselves with a certain degree of model dependence 
in the data. Realistic tests of the Copernican principle will 
therefore have to take special care to properly model the po¬ 
tential function $ employed to perform inference. Another 
possibility is to use the model dependent data sets only as a 
means to specify priors on certain functions and then use the 
model independent data to perform inference. These consid¬ 
erations will be addressed in future research. 

4.3 Implementation 

In order to write down the final form of the algorithm we 
need to clearly state our objectives. The primary aim of this 
work is to test the Copernican principle. In order to do so 
we need to find the full distribution of solutions that are 
compatible with currently available data and then use these 
to reconstruct the quantities Ti and T 2 defined in 2.2. How¬ 
ever we still need to decide in advance on the domain over 
which the solutions should be found. Since the domain in v 
is completely determined by the observations this amounts 
to specifying a time in the past up to which to integrate to. 
However this choice is not straightforward for two reasons. 
Eirstly, the time is not arbitrary since, as explained in sec¬ 
tion 2.3, W will intersect C at a finite time in the past. The 


available data therefore also places a limit on how far back 
in time we are able to go. Secondly, since there is also signif¬ 
icant variation in the current age of the universe to, we do 
not simply want to choose a “safe” value that will always lie 
within the causal horizon. The reason for this is that it might 
happen that a large percentage of the models considered end 
up having a to smaller than the value we choose to integrate 
to. We therefore need to find a compromise between these 
two extremes. Eor this paper we have chosen to integrate up 
to t* ~ 10.5 Gyr. Henceforth we refer to the PLC defined 
by = t* as PLCE. We will also be reconstructing the dis¬ 
tributions of JA, R and p on the hypersurface St defined by 
t* — const. 

Having specified the domain we can find the full distribution 
of solutions inside the PLC up to the PLCE using Algorithm 
1. Eor this algorithm we assume that initial realisations of 
H{z)^ p{z) and A have been chosen such that they satisfy all 
the physical requirements of section 3.2. The potential func¬ 
tion is then initialised by computing (53) with the resulting 
model. 

The code implementing Algorithm 1 is available on request 
and will be made publicly available at a later stage. It has 
been written in a mixture of Python and Eortran. This code 
was run on the Rhodes maths cluster which is hereby dubbed 
the GetaFix cluster. 

4.4 Simulations and verification 

To test our numerical implementation of the above algorithm 
we simulate data around a background ACDM model defined 
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Figure 6. Observables on the PLCO. Top Left: Posterior distribution of D(z) on the PLCO for real data out to 2 ; = 1.75. Top Right: 
Posterior distribution of H{z) on the PLCO for real data. Bottom Left: The posterior distribution of p{z) on the PLCO for real data. 
Bottom Right: The posterior distribution of to- The black (green) cross indicates the background ACDM (LTB) value. The large degree 
of uncertainty in to mainly stems from the uncertainty in A. This results in a small percentage of the initial samples having a value of 
to < t*. For comparison we show, in green, these relations for a typical LTB model. 


Algorithm 1 

1 Propose + /3(5, (5 ~ /io(ic) and 

sample A ; 

2 if any(p < 0^ then $ = 00 go to 6; 
else Implement Procedure 1; 

3 if to < t* then $ = 00 go to 6; 
else Implement Procedure 2; 

4 if >V intersects C at w > t* then $ = 00 go 
to 6; 

else Implement Procedure 3; 

5 if Shell crossing then $ = 00 go to 6; 
else Compute ^(^(fc)) with (53); 

6 Set X(^k+i) — y{k) with probability (52) ; 

7 Set k ^ k 1 and go to 1; 


by 

^mo = 0.3, Qao = 0.7, iLo = 70 km s“^Mpc“^. (55) 

To simulate the data we have followed a similar procedure to 
that described in Bester et al. (2014) except that we chose 
a simple power law relationship for how relative error scales 
with redshift viz. 

6{z)^6o{l^zy. (56) 

Here (5o is the relative error of the simulated data sets at 
z = 0. For both D and H we chose | with a relative er¬ 
ror at the origin of 5% and 10% respectively. For the mock 
density data we chose a = 4 with a relative error at the 
origin of 40% and error bars centred on p. We simulate a 


total of 500 data points for D and 50 for both H and p with 
redshift values drawn from a uniform distribution between 
z — 0 and z = 2.0. Finally, as explained in section 4.1, we 
use a Gaussian prior for A with mean corresponding to the 
background ACDM value and an uncertainty of 2%. 

We then run these data through Algorithm 1. To make sure 
that the distributions have converged we run a number of 
these processes in parallel, each with a 10% burn in period. 
The first process draws 1000 samples and each subsequent 
process doubles the number of samples generated by the 
last. We keep spawning new processes until the posterior 
distributions stops changing. Typical acceptance rates vary 
between 25% and 35%. 

The results are summarised in Figures 3-5. In all these fig¬ 
ures the black crosses indicate the background ACDM val¬ 
ues. Contours are found by constructing the empirical dis¬ 
tribution function at each point in the domain. In all these 
figures the blue line is the median, the lightblue (dark blue) 
shaded area is the 1-a (2-cr) contours as they contain 68% 
(95%) of all realisations. As can clearly be seen from the fig¬ 
ures the background model always falls within at least the 
2-a contours of the reconstructed functions but are more of¬ 
ten confined to the l-cr contours. It is especially reassuring 
that neither Ti nor T 2 show a statistically significant depar¬ 
ture from zero in Figure 4. Also note that the most probable 
density profile in Figure 5 is very nearly flat but that void 
profiles do not seem to be ruled out. Figure 5 also shows that 
the metric functions X and R, although they exhibit signif¬ 
icant uncertainty, both agree with what is expected from a 
flat ACDM model i.e. that X is constant and R a linear 
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Figure 7. Tests of the Copernican principle. Top Left: Posterior distribution of Ti on the PLCO marginalised over to- Middle Left: 
Posterior distribution of Ti on the PLCF. Bottom Left: The maximum extent of v on the PLCF. The black cross is the background 
ACDM value. Top Right: Posterior distribution of T 2 on the PLCO marginalised over to- Middle Right: Posterior distribution of T 2 on 
the PLCF. Bottom Right: The maximum extent of 2 ; on the PLCF. The black (green) crosses are the background ACDM (LTB) values. 
Clearly ACDM is compatible with the data. However we cannot yet rule out LTB model shown in green. 


function of the comoving distance r. 

There seems to be a very good agreement between the recon¬ 
structed functions and the background ACDM values. This 
confirms that our numerical implementation of Algorithm 1 
performs as expected. 


5 RESULTS AND DISCUSSION 

The result of running currently available data through Algo¬ 
rithm 1 are summarised in Figures 6-8. For comparison we 
include a typical LTB model in the figures. The LTB model 
we use is the constrained LTB model of Garcia-Bellido & 
Haugboelle (2008) with parameters chosen to closely mimic 
the ACDM D{z) relation. We are able to go out to a red- 
shift of 2 : = 1.75 as that is the maximum redshift for which 
we have H(z) data. This places constraints on D(z) be¬ 
yond the maximum redshift of z = 1.41 in the Union 2.1 
sample. Note that we use the same prior over p(z) as for 
the simulated data. As explained in section 4.1 we use a 
flat prior for A. These results confirm that the data are 
compatible with ACDM. Moreover we find the value of 
Hq = 69.76io!7i kms“^Mpc“^ at the 2-a confidence level. 
However the tight constraints on Hq are artificial and result 
from not marginalising over its value when converting dis¬ 
tance modulus to angular diameter distance data. 

The similarity of the reconstructed distributions of T 2 on the 
initial and final PLC’s, for both simulated and real data, sug¬ 
gests that it does not change by much during the evolution 
of the universe, at least for the class of models considered 
in this work. Furthermore this quantity seems to be quite 
insensitive to the value of A. It may therefore be able to 


shed some light on the Copernican principle even with our 
current uncertainty about this elusive parameter. It should 
be noted that the constraining power of angular diameter 
distance and expansion rate data diminishes at high red- 
shifts. Density data at high z should significantly constrain 
T 2 . However, given that realistic models of the universe (i.e. 
FLRW + perturbations) should have a value of T 2 as low as 
1 ^ 2 ! 10 ^ Maartens (2011), it remains to be seen whether 

density data from upcoming surveys will be sufficiently ac¬ 
curate to test the Copernican principle. A more detailed 
investigation is in preparation and will be released at a later 
stage. 

The value of A seems to have a greater impact on the evo¬ 
lution of the quantity Ti. However the reduced uncertainty 
of Ti as compared to T 2 at high redshifts suggests that Ti 
might be more appropriate to test the Copernican principle 
on large scales. Currently available data does not place very 
stringent constraints on the relative difference between H± 
and H. On the PLCO these two quantities remain, at the 2-a 
confidence level, within 10% of each other up to z ~ 0.15 
but can differ by as much as 30%-40% ai z ^ 1.75. As can 
be seen from Figure 7 these constraints quickly degrade on 
the PLCF. Figure 4 shows that the constraints are signifi¬ 
cantly better for our simulated data. Since the uncertainties 
in D{z),H{z) and p{z) are very similar for the simulated 
and real data it is the uncertainty in A that degrades the 
constraints on Ti. Redshift drift data should therefore con¬ 
tribute significantly to tightening the constraints on Ti. This 
is another avenue that we are currently investigating, further 
details will be released at a later stage. 

Neither Ti nor T 2 are able to reject the void LTB model used 
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Figure 8. Posterior distributions on a constant time slice St at = 10.5 Gyr for real data. Top Left: The metric function r*). This 

function is compatible with what is expected in ACDM since it is very nearly constant. However, because of the small volume probed by 
the current observations, we see that the LTB model (shown in green) seems to follow the same trend. The data would have to probe a 
much larger volume for the departure from this trend to be visible to the naked eye. Bottom Left: The density on St plotted in units 
of the critical density pc of the background ACDM model used for the simulated data. Once again the most probable density profile is 
flat. However void profiles are not ruled out by the data. The density profile of the LTB model shown for comparison correspondes to 
that of a very shallow void. The overall value of this density is significantly smaller than that of ACDM and is almost ruled out at the 
2-(j confidence level. Top Right: The metric function R(t*,r*) is very nearly a linear function of the comoving distance r and is thus 
compatible with FLRW. However the same remarks as those regarding X* apply here as well. Bottom Right: The maximum extent of 
the radial coordinate r* on St. The black (green) cross is again the background ACDM (LTB) value. 


for comparison. In fact the LTB model lies comfortably in¬ 
side the 1-cr contours of the current data. The observational 
signature of this model on the quantity T 2 is especially small. 
We found that using H = H± in (20) results in a larger de¬ 
viation of T 2 from zero. However the constraints from the 
data similarly degrade under this substitution so we do not 
gain anything from it. The inability of these tests to reject 
the LTB model suggests that much better data are required 
to test the Copernican principle. 

Finally in Figure 8 we show the distributions of the met¬ 
ric functions X and R, as well as the distribution of p, on 
Et. We have again included the LTB model for comparison. 
Both metric functions seem to follow trends typical of a flat 
FLRW cosmology i.e. that X is constant and R is a linear 
function of r*. Although there is significant uncertainty in 
the values of these quantities it is not their values that we are 
after. We are really more interested to find out whether these 
functions can be considered as separable functions of t and 
r. If they can it would be a nail in the coffin for radially inho¬ 
mogeneous models since any spherically symmetric universe 
with separable metric components has to in fact be maxi¬ 
mally symmetric and therefore homogeneous. Unfortunately 
it is very difficult to say with statistical certainty whether 
this is indeed the case. In fact the LTB metric components 
shown in these figures seem to follow the same trends. This 


is because of the relatively small volume that is accessible 
through observations at t = 10.5 Gyr. From the figure we 
see that this volume is less than 1 Gpc^ whereas the full 
width at half maximum height of the void used to generate 
the LTB model is approximately 1.5 Gpc at t = to- The near 
flatness of the LTB density profile shows that the void has 
be very shallow at early times. However we also see that the 
LTB model’s density would have been considerably smaller 
in the past so it might be possible to rule out certain classes 
of LTB models based on the allowed values of the density. 
Unfortunately we cannot draw any strong conclusions from 
this fact alone. 

We have presented an algorithm that can in principle con¬ 
front ALTB models with arbitrary data. Moreover the func¬ 
tional degrees of freedom of the ALTB model are fixed us¬ 
ing directly observable quantities that are free from any 
gauge effects. We have further illustrated how A, the one 
free parameter of the model, can be constrained directly 
from observations of the redshift drift which is in principle 
completely model independent. As such we have presented 
a general framework for testing the Gopernican principle 
both with currently available data and data from upcoming 
surveys. We have intentionally refrained from incorporat¬ 
ing GMB data. The main reason for this is that the dy¬ 
namics of non-comoving dust and radiation fluids in spher- 
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ical symmetry is poorly understood Lim et al. (2013). Self 
consistently incorporating radiation into the observational 
cosmology programme is highly non-trivial because we do 
not want to presuppose that the early universe is homo¬ 
geneous. Instead the ideal is to model late time processes 
as accurately as possible and constrain the geometry of the 
universe based on data gathered with these models. Such 
an approach would be complementary to that employed in 
Redlich et al. (2014) and Zibin & Moss (2014) for example. 
Both approaches should ultimately arrive at the same con¬ 
clusions. 

Furthermore in section 4.2 we suggested a way to incorpo¬ 
rate data that are not entirely model independent. This is 
obviously not the ideal way to test the Copernican principle. 
However the considerable difficulty involved in getting model 
independent data might necessitate this kind of approach for 
some time to come. We should also keep in mind that this 
framework relies on the validity of a number of additional 
assumptions. As a result the tests we perform do not only 
test the Copernican principle but also the compatibility of 
the data with these assumptions. We should therefore strive 
to confront the data with as many consistency relations as 
possible. It would seem, at least within the framework pre¬ 
sented in this paper, that there is some work to be done 
before we can ultimately conhrm or refute the validity of 
the Copernican principle. 
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